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Abstract 

We compute the production of gluons from Glasma color flux tubes. We calculate the probability 
distribution of gluon multiplicities arising from the distribution of color electric and color magnetic 
flux tubes found in the Glasma. We show that the result corresponds to the negative binomial 
probability distribution observed in experiments. The parameter k that characterizes this distri- 
bution is proportional to the number of colors N c 2 — 1 and to the number of flux tubes. For one 
gluon color and one flux tube, the multiplicity distribution is close to a Bose-Einstein distribution. 
We call this decay process "Glitter" , a term that is explained below. 

Key words: glasma, multiplicity distribution 
PACS: 13.85.Hd, 24.85.+p, 25.75.-q 

1. Introduction 

In high energy nuclear collisions, it has been argued that a Glasma is formed through the 
collision of two sheets of Color Glass Condensate (CGC) Q, 0, 0, 0, 0, 0, E3, El El, E3, 

14, 15, El; 12, 03] (for a review and additional references see [19L [20f| ) . At very high energies, 
the gluon density per unit area per unit rapidity, dA r /dyd 2 rj_, is of order Ql/a s in the CGC and 
Glasma. This density is so large that the interaction strength of QCD is weak, a s < 1. Indeed, 
in these relations, the scale of the coupling constant is set by the saturation momentum Q s , and 
the saturation momentum grows with both increasing energy and size of the nucleus. The Glasma 
is formed during the time it takes two Lorentz contracted sheets of Colored Glass to pass through 
one another, r ~ e~~ K / as /Q s a time parametrically short compared to the natural time scale for 
decay of the flux tubes, Td eC ay ~ 1/Qs- We shall not describe in detail the properties of either the 
Glasma or the CGC in this paper and refer the reader to the original literature for details. 

In this paper, we compute the probability distribution for the multiplicity of gluons produced 
in the Glasma. We show that the distribution of gluons arising from such a decay is not a Poisson 
distribution as might be expected for the decay of an external source into particles in a weakly 
interacting theory. It turns out the distribution is a negative binomial distribution. 

Recall that a Poisson distribution, 

1 
i\ 

is completely characterized by its mean value n. In contrast, a negative binomial distribution is a 
2-parameter distribution of the form 

nb _ T(k + n) n n k k 

n ~ r(fc)r(n + l) {n + k) n + k ' ( ' 
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pPoisson = -^n n e- n (1) 



This distribution has larger fluctuations than a Poisson distribution, but tends to a Poisson distri- 
bution if k — ► +00 at fixed n. 

For the processes we consider, the decays will be into one coherent state associated with a 
gluon. For this reason, and because a negative binomial distribution does not fall as 1/n! at large 
n, as does a Poisson distribution, we will refer to this property of the distribution as tenacious. The 
tenaciousness of this distribution results in an amplification of the intensity of multiply emitted 
gluons relative to that of a Poisson distribution. This means that there are larger fluctuations in 
high gluon multiplicity events than would be typical of a Poisson distribution. Thus, we shall use 
the acronym "Glitter" to describe radiation from these flux tubes, as an abbreviation for GLuon 
Intensification Through Tenacious Emission of Radiation. 

In this paper, we shall compute the multiplicity distribution of gluons produced by Glasma flux 
tubes. We show that a single flux tube decays into gluons with a negative binomial distribution 
and is characterized by a parameter ko of order one. This implies that a decay of N FT flux tubes 
produces a negative binomial distribution characterized by k = N FT ko. For k = 1, a negative 
binomial distribution is a Bose-Einstein distribution^, so single flux tube decays are close in form 
to Bose-Einstein distributions. In addition to the parameter k, the negative binomial distribution 
is also parameterized by the average multiplicity n. The ratio n/k = n/N FT ko is approximately 
the multiplicity per flux tube. 

Based on this interpretation of the decay of an ensemble of flux tubes we argue that the flux 
tubes are a "glittering" glasma. We will then iscuss our results in the context of the extraction of 
negative binomial distribution from the UA(5) and PHENIX collaborations. 

2. Calculation of the multiplicity distribution 

The probability distribution of a discrete quantity is conveniently defined in terms of its gen- 
erating function 

00 

F(*) = J>"P n , (3) 

n=0 

from which one can compute the moments of the multiplicity distribution as 



(n(n-l)---(n-q+l)) - d<?F(z) 



dzi 



(4) 

z=l 



(The moments defined in this way are known as the factorial moments.) It was shown in Refs. 21 



22j, |23j] that in the case of the central rapidity region of nucleus-nucleus collisions, when both nuclei 



can be described as strong color sources p ~ 1/g, these moments can be computed as 

( n(n _ 1) . . . ( n _ q + i)) = j [d Pl }[dp 2 ]W y [ Pl ]W y [p 2 ] {n[p 1 ,p 2 ]) q . (5) 



This result is valid to leading log accuracy, i.e. it includes the leading order in a s with all the 
powers of a s In 1/x resummed into the rapidity dependence of the weight functionals VFyl/Oi]. The 



X A Bose-Einstein, or geometrical, distribution is a thermal distribution for single state systems. Its probability 
distribution reads: 

be 1 / n 

1 + n U + n , 
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Figure 1: The fundamental building block, the one gluon production amplitude. 
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Figure 2: Diagrams that have to be contracted into connected one. 



factor n[p\, P2] inside the integral in the r.h.s. of eq. © is the integrated multiplicity corresponding 
to a fixed configuration of color charge densities p\ and p2- 



n\Pi,P2\ 



d P±dy p - 



dN 



(6) 



d 2 p±dy p ' 

where dN/d 2 p^dy p is obtained by Fourier transforming the classical gauge field radiated by the 
sources p\p- In our leading log calculation the classical fields and thus also the single gluon 
multiplicity are boost invariant, so the integral over rapidity in eq. ([6]) is just a constant factor. 
We are assuming that the rapidity interval is small enough compared to l/a s ; otherwise there are 
additional large logarithms that must be resummed; this case is studied in detail in Ref. j23| . 

In the Glasma the single inclusive multiplicity is of order 1 /a s and thus the moment defined in 
eq. © is of order (l/a s ) q . We are computing the moments only to leading order in a B , thus powers 
of n lower than q are negligible compared to the g'th moment. In particular (n(n— 1) • • • (n—q+1)) ~ 
(n q ) at this level of accuracy. This means that the leading (in a s ) correlations in the multiplicity 
distribution come entirely from the average over the distribution of sources pi 2) i- e - the dominant 
behavior of the probability distribution comes from the large logarithms of the energy resummed 
into the W's. The contributions to the probability distribution for fixed sources from higher loop 
orders that were studied in Ref. [24}] are suppressed by powers of a s and thus contribute to our 
calculation only when they are enhanced by large logarithms. 

The fundamental properties of the probability distribution are better reflected in the factorial 
cumulanta^l 

d 9 lnF(z) 



m q = (n(n — 1) • ■ ■ (n — q + 1)) — disc. 



dzi 



(7) 



where "disc." denotes the disconnected contributions that can be expressed in terms of the lower 
cumulants. We shall now turn to calculating these quantities of the multiplicity distribution in the 



2 Note that the factorial cumulant defined here differs from the conventional one in that we are defining the 
cumulant from (n(n — 1) • • • (n — q + 1)) instead of (n q ). In terms of the generating function this is a question of 
differentiating w.r.t. z in stead of \nz. As explained above the difference between the two is of higher order in a s 
than our calculation and we will not discuss it further. 
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Figure 3: Contraction contributing to the dominant correlation, building block of rainbow diagram 




Figure 4: Contraction contributing to a subdominant correlation, non-rainbow diagram 



Glasma at the lowest nontrivial order in the sources. This approximation is equivalent to assuming 
that the momenta of the produced gluons are all larger than the saturation scale. For concreteness 
one can take the distribution of the sources from the MV model: 



W[p)=C exp 



2 p a {x_ L )p a (x_ L ) 

d x± t-= 



(8) 



Since the probability distribution essentially depends only on the combinatorics of pairwise source 
connections, our result applies equally well to a nonlocal Gaussian distribution that would more 
closely reproduce a solution of the BK equation 25, 0,0- 

The calculation proceeds in the same way as that of the second and third cumulants computed 
in Refs. 



271 . 1281 . l29(, and we refer to these works for a more detailed description. Computing the 



probability distribution to all orders in the sources /? i p, is in principle possible using the methods 
developed for the single inclusive gluon production pj], EH, Efl, 17], but reliably computing the 



higher cumulants requires a significant numerical effort to gather enough statistics and is left for 
future work. 

The fundamental building block in our calculation is the amplitude to produce one gluon with 
momentum p from the fixed classical color charges pi(k±) and P2(p± — k±). This amplitude reads 



Pl{k±) p2(p± ~ 
{p ± -k 



k\ 



(9) 



where we are not writing the color indices explicitly. Here, -L 7 (p, denotes the effective Lipatov 
vertex. We do not need the explicit expression of its components, and it will be sufficient to know 
that it satisfies the following two properties 



r<( P ,k A 

L 7 (p,fc_L)L 7 (p,fcj 



£ 7 (p,P_l - k±) , 
, ( P± -k ± ) 2 kl 

P± 2 



(10) 



Figure 5: Rainbow diagram. 
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Figure 6: Dimer notation for rainbow-like links. 



The diagrammatic notation for this amplitude is shown in fig. [TJ To compute the g'th cumulant 
we need to take 2q factors of this basic building block (q for the amplitude and q for the complex 
conjugate) and perform the averages over the sources. Because the distribution of the sources in 
eq. ([H]) is Gaussian, we only need to keep track of contractions of pairs of sources p\ and (separately) 
of pairs of sources p2 and replace them by the correlator 

(p(k ± )p(k' ± )) = (27r) 2 5\k ± + k' i _)g i p 2 {k 1 _). (11) 

For the time being we shall leave an unspecified fc^-dependence in the correlation function g 4 p 2 (k±). 
In the MV model [H, H, 0] g A u 2 (k \) i s a constant, but JIMWLK or BK evolution can effectively 
lead to a different fc^-dependence [301 ]. With the simplified diagrammatic notation introduced in 
fig. [T] this combinatoric problem now corresponds to forming a connected contraction of the 2q 
boxes, each with two lines attached, illustrated in fig. [2j 

We shall now show why the only contributing contractions are "rainbow" diagrams, where on 
one side (upper or lower) of the diagram the two boxes corresponding to the same momentum p r 
are contracted with each other, as in fig. [SJ An example of a building block of a rainbow diagram 
is given in fig. [3l There are a total of four propagators with the same momentum k±. Two of 
these are cancelled by the contractions of the four Lipatov vertices attached to the ends, leaving 
a quadratically infrared divergent contribution ~ d?k^_/k\_ to the integral over k±. These kinds 
of divergences are a sign that this contribution to the multiparticle correlation is sensitive to the 
whole correlated area in the transverse plane. They are regulated at the scale Q s (since Q~ x is the 
correlation length between color charges in the transverse direction), giving a contribution of order 
1/Qs- Compare this to the "non-rainbow" contribution depicted in fig. HI Here there are only two 
propagators with the same momenta and the fcj_-integral is convergent. Instead of 0(1/Q 2 ), this 
yields a contribution ~ 1/P±i or ~ 1/Pj_2 which shall be neglected here since we are assuming 
p_i_ 3> Q s . Thus the contractions of the boxes in fig. [2] have to form a "rainbow diagram" (fig. [5]) 
in either the upper or lower part of the diagram. There cannot be a rainbow on both sides since 
this would lead to a disconnected contribution. 

Now that we have reduced the combinatoric problem to rainbow diagrams we can simplify our 
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Figure 7: Connected diagram in the polymer notation introduced in fig. [6] 



diagrammatic notation even further, as shown in fig. El We can consider the two boxes connected by 
the line in the rainbow as a dimer with two ends. Now we must count the number of ways in which 
our q dimers can be combined to form a connected loop. The first end of the first dimer can be 
connected to 2q — 2 other loose dimer ends, let us assume it is connected to dimer r. The second end 
of dimer r now has 2q — 4 loose ends available, because closing the loop with the first dimer would 

immediately give a disconnected contribution. Thus we arrive at (2g-2)-(2g-4) 2 = 2 q - 1 {q-l)\ 

combinations. Remembering that we can form the rainbow on either the upper or lower side of the 
diagram we finally arrive at 2 q (q — 1)1 topologies. In the special cases of q = 2 and q = 3 this gives 
4 and 16, which were the number of diagrams evaluated in 27] and [2^ ]. 

The color structure can be simplified in the following way. Let us denote the color indices of 
the sources i = 1 . . . q contracted on the upper side of the diagram in fig. [5] by 04 (because of the 
rainbow structure of the connections on the upper side, we do not need separate indices for the 
sources in the complex conjugate amplitude), the color index of the gluon with p i by bi and the 
color indices on the lower side of the diagram Cj (source connected to gluon p i in the amplitude) 
and (in the complex conjugate). Now the color structure of the upper side of the rainbow and the 
boxes is fo-ibm faibic[ _ _ _ ja q b qCq p q b q c' q _ (c A )v8 c i c 'i . . . S c " c 'i. C A = N c is the Casimir of the adjoint 
representation. The Cj,c^ indices must now be contracted pairwise into a single connected loop as 
in fig. which yields a factor tr(l a( jj) = iV c 2 — 1, making the total color factor (N c ) q (N c 2 — 1). 

Let us then turn to the structure of momentum flow in the diagram. Transverse momentum 
is conserved at the vertices and in the sources connections due to the expectation value eq. (fTTj) . 
Altogether there are originally Aq transverse momentum integrals from the powers of the sources. 
There are 2q delta functions from the source correlators eq. (jlip and 2q momentum conservation 
delta functions from the three gluon vertices. Not all these delta functions are independent: two 
of them end up having the same argument, yielding one factor of the transverse area (denoted 
S±). Therefore, there is only one remaining transverse momentum integral. One can choose this 
remaining momentum to be the one circulating in all of the lower part of the diagram, which we 
shall denote by k± (this is the momentum that circulates along the loop in fig. E]). On the rainbow 
side of the diagram there is a squared propagator l/(k± — p± i ) 4 for all the sources i = 1 . . . q. On 
the non-rainbow side the transverse momentum in all the propagators is the same, giving a factor 
l/k%. Half of these propagators are cancelled by the squares of the Lipatov vertices, which also 
contribute an inverse square of the external momentum. Combining the combinatorial factors from 
the source averages, the propagators, Lipatov vertices and factors from the invariant measure, we 
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finally getH 



dN 

dyid 2 p ±1 . ..dy q d 2 pj 



2 q (q-l)l 



d 2 k ± fg 4 ^ 2 (k 
(2vr) s 



(N C )"(N C 
(P±i) 2 ■ ■ ■ (P± q 



1)5, 1 



2<i 



j_ 

I 2 g 2 i{2n)^ 



k l 



(Pj.i 

IE 



(Pi 



A; 



(12) 



for the single inclusive spectrum 



This general formula also reproduces the result of refs. 
case q = 1; in this case the combinatorial factor in the square bracket must be taken to be 1 instead 
of 2 to avoid double counting the only contributing diagram. 

The weak source result eq. (|12|) is infrared divergent in the MV model (g 4 fJ> 2 (k±) constant). 
Physically this is modified by several effects. Even in the weak field limit BK or BFKL evolution 
leads to an anomalous dimension < 7 < 1 that changes the behavior into g 4 fj. 2 (k±) ~ fc^ 1 7 ^ 
in the geometric scaling region k± > Q s . Deep in the saturation regime it has been argued 
that the correlator effectively behaves as g A [i 2 (k±) ~ k\. Ultimately the infrared behavior of the 
multigluon spectrum is regulated by the nonlinear interactions that are not included in our present 
computation. This is seen explicitly and analytically in the "pA" case [33, 34, 35] and in numerical 
computations of the glasma fields in the fully nonlinear case 14, 15, 1|| 17]. Since the full nonlinear 
dynamics are known to regulate the infrared behavior in the case of the single gluon spectrum we 
have strong reasons to expect that they will also do so in the case of multiple gluon production; 
at the same scale k± < Q s . We emphasize that an essential point in this argument is that the 
quantity appearing in eq. (|12p is not a single color charge correlator divided by a large power k 2 ^, 
but the same correlator g 4 ^ 2 (k±)/k 2 L that appears in the single inclusive gluon spectrum raised to 
a large power q. 

The effect of saturation on the multigluon spectrum at k± < Q s has a very intuitive interpreta- 
tion in the the glasma flux tube picture. The size of the flux tube, 1/Q S , is the correlation length of 
the system and we should not have contributions from longer distance scales. We effectively take 
this into account by regulating all the infrared divergences at the scale Q s , and thus approximating 

the integral in eq. (|12|) by (^2tt Q 2 ^ q 2 ^p± 2 ■ ■ ~P± 2 ^j ■ Here k is a constant of order one that 
depends on the details of how the infrared divergences are regulated at the scale Q s . Since Q s is 
the typical momentum of the produced gluons, not a lower limit, we expect that numerically k < 1. 
In our analytical calculation we do not have access to the exact value of this coefficient. We also 
use the corresponding approximation for the single inclusive spectruixE]: 



diV 
dyd 2 pj 



N C (N C 2 - 1) S L (g 2 f,) 4 



4vrV 



P± 



(13) 



3 In comparing this to the result in 27| (eq. (27)) note that there is an erroneous factor of 1/2 in eq. (9) of [2t| 
which propagates into an additional factor 2~ 2q . In addition, there is an overall error of (2tt) 2 . These errors were 
present also in Ref. but have been corrected in the published version. 

4 We are neglecting an additional logarithmic dependence in p ± . 
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We can now express our result as 



diV \ _ (N c 2 -1)kQ 2 S ± ({g 2 ^ 1 N c \ 9 1 



d yi d 2 P±1 ...dy q d 2 p ± J v " r 2vr V 9 2 I^kQI) (p ±1 )4 • • • (p, )4 



9" 

dJV \ / dN 



W j ' 2vr ((iV c 2 -l) K Q2 5±/(27r)) « I ) 

If we integrate this equation over the rapidities and transverse momenta of the q gluons, again 
consistently regulating all the infrared divergences at the scale Q s , we obtain our result for the 
factorial cumulant as 

'n\i 



m q = ( q -l)\k{-) , (15) 
with 

The exact constant factors, encoded in the coefficient k, depend on the exact way the infrared 
divergences (logarithmic for the single inclusive, power law for the multigluon correlations) are 
regulated. These factors cannot be obtained exactly in an analytic calculation to the lowest order 
in the sources. However, the main parametric dependences in the relevant variables a s ,Q s ,S± 
and N c can be expected to be the same to all orders in the sources. A possible additional (mild) 
(/-dependence in n would be a minor correction to the behavior of the probability distribution, 
mostly determined by the combinatorial factor (q — 1)!. 

Equations (|15p and (|16p are the main result of this paper. One can see that these factorial 
cumulants (|15p are those that define the negative binomial distribution. It arises very naturally 
in the Glasma based on the Gaussian combinatorics of the classical sources and the assumption of 
the fluctuations in the system being dominated by a correlation length 1/Q S . 



3. Glittering Glasma: Interpretation of the result 

A negative binomial distribution is characterized by two parameters, the mean h and k, in 
terms of which the probability to produce n particles is 



NB 



Y{k + n) n n k 



Pn T(k)T(n + 1) (n + k) n + k ' ^ 



The distribution is characterized by the generating function 



F ktfi (z)=Y / z n P n =fl-l(z-l)) . (18) 



k 

n=0 



The moments of the distribution can be obtained from the generating function by differentiating 
with respect to z. The connected parts of the moments, or cumulants, are generated by the loga- 
rithm of the generating function. The factorial cumulants ([7]) of the negative binomial distribution 
are given by 



d q 

- — In Fu (z) 

dzi fe ' nV ; 



(g-l)Ifc(?V. (19) 



In contrast, for a Poisson distribution mi = n and m q = for q > 1. This quantity is the 
expectation value 

m q = (n(n — 1) • • • (n — q + 1)} — disc, (20) 

where the disconnected part "disc" can be expressed in terms of the lower order cumulants. The 
convention that the expectation value in ([20]) is taken of the product n(n — 1) • • • (n — p + 1) and 
not of n p means that we are subtracting the "Poissonian" part from the moment, which is why 
all the factorial cumulants of a Poisson distribution are zero for q > 2. Note that the Poissonian 
part is suppressed by powers of a s in the Glasma; thus in practice the difference between m p and 
a conventional cumulant is neglected in our analysis. 

Two common special cases of the negative binomial are the Poisson distribution, obtained in 
the limit k —* oo at fixed n, and the geometrical-or Bose-Einstein-distribution obtained when 
k = 1. The negative binomial distribution is wider than a Poisson distribution typically associated 
with independent emission of particles; this can be seen e.g. from the variance 

a 2 = (n 2 )-(n) 2 =n + ^-. (21) 

A useful property of the negative binomial distribution with parameters n, k is that it is is also 
the distribution of a sum of k independent random variables drawn from a Bose-Einstein distribu- 
tion with mean n/k. This is easily seen from the generating function in eq. (I18p . remembering that 
the generating function of a sum of independent random variables is the product of their generat- 
ing functionqf|. This has a consequence that an incoherent superposition of N emitters that have 
a negative binomial distribution with parameters fco , n produces a negative binomial distribution 
with parameters Nko,Nn. 

A natural physical interpretation of our result can be given in terms of emission from indepen- 
dent glasma flux tubes. Geometrically, the transverse area S± is filled with Q 2 S± independent flux 
tubes of size ~ 1/Qs- Each of these tubes emits gluons in N c 2 — 1 different colors. Our result shows 
that the probability distribution of gluons of one color emitted from one flux tube is approximately 
a Bose Einstein distribution. We do not see an interpretation of this result as a thermal process, 
however. It seems more likely that the distribution is one that maximizes the entropy (which is 
the defining property of the BE distribution) because there is a large number of color sources that 
emit gluons. In some sense the role of a heat bath (large reservoir of energy) in thermodynamics is 
played by the large number of color charges resummed into the effective color current of the CGC. 

It is a known experimental observation 3^, 3?], 2^ (see also 40] for an extensive review) that 



multiplicities of charged particles in high energy scattering are well described as a negative binomial 
distribution. Also multiplicity fluctuations at RHIC have been found to agree with the negative 



binomial distribution by the PHENIX collaboration [4l|, |42|]. Experimentally, the parameter k 
increases somewhat with 5r], the size of the rapidity interval in which the particles are measured. 
The dependence is, however, very slow for large 5r], pointing to the presence of a long range 



correlation in the system [43]. This is natural in the Glasma picture, since flux tubes extend over 



large rapidity intervals. The number of flux tubes, which gives the parameter k of the negative 



5 Explicitly, consider n = ni + ■ ■ ■ + n r where the rii's are independent of each other. The probability dis- 
tribution of n is then P n = Y2 ni '"53 n ~ n 0-P«i ' ' ' Pn r an d the generating function ^2„z n P n — 
X] ni ' ' ' X^n hnr ^n 1 ■■■Pn r , which is the product of the individual generating functions for the variables 
rii. 
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45] 
411 



binomial distribution, essentially depends only on the transverse area of the projectiles and on the 
saturation momentum. 

The main difficulty in interpreting the experimental results arises from the geometrical fluctu- 
ations from averaging over different impact parameters in one finite centrality bin. To minimize 
this effect one should use as small centrality bins as possible. Comparison with a different method 
of analysis used by the STAR collaboration 44| could be very useful in disentangling these ef- 
fects. To the extent that this uncertainty allows us to compare results in gold-gold and pp, the 
picture we present seems fairly consistent. For a fixed collision energy we would expect sca ling 
k ~ QgS± ~ iVpart- While keeping this caveat in mind, the results from UA5 [1^] and E735 
(Jfe « 2 ... 4, A/p ar t = 2) and PHENIX k ^ 350 for 0-5% most central ( A part « 350) collisions 
or k = 690 when extrapolated to a zero centrality bin width 42|, seem very consistent with this 
estimate. 

From the experimental fit of the parameter k in central gold-gold collisions by PHENIX [13] 
and the value Q s ~ 1.1 GeV estimated from measurements of the charged multiplicity, one can 
use eq. (|16p to obtain an estimate n ~ 0.2 for the parameter that reflects our uncertainty in the 
infrared sector. As we have discussed, it is natural to expect a numerical value of slightly less 
than 1 for k. This, however, means that at RHIC energies the flux tube size, as measured in the 
multiplicity distribution, is not yet very clearly separated from the confinement scale. At LHC 
energies we can expect this separation to be clearer. 

For increasing collision energy we would expect Q s and therefore k to increase. The energies 
where the UA5 measurements are done are still in the transition region from a behavior of k 
decreasing with energy from lower y/s, but we would expect k at the LHC to be clearly larger. 
This decreasing behavior at low energy follows because of the Poisson nature of low energy particle 
emission, and that for a Poisson distribution k — > oo. 

The negative binomial has been interpreted as resulting from a partial stimulated emission 



or cascade process [43J . It has been known in the literature [46j, |47|] that the distribution would 
naturally arise from a superposition of subsystems with Bose-Einstein distributions. Nevertheless, 
a popular app roach has remained to interpret the observations in terms of a fluctuating number 
of strings [43], each producing particles typically with a Poisson distribution [i^, (see also 
l53| ] for a more pQCD based approach). While the picture of flux tubes in the glasma has 



many similarities to ideas in string model phenomenology, the distribution of particles produced 
from one flux tube is different. The probability distribution of gluons from a glasma flux tube is 
not a narrow Poissonian, but has very large fluctuations: the glittering of the glasma. 



4. Summary 

The Glasma provides a successful phenomenology of a particle production in high energy 
hadronic collisions. There is now experimental data on the ridge phenomena that show flux 
tube structures in two particle correlations 54, I55I IslUl . In addition, long range correlations of 
remarkable strength are seen in heavy ion collisions [441 ] . 

The Glitter of the flux tube decay may provide a strong tool for disentangling various de- 
scriptions of the flux tubes, since it naturally leads to a negative binomial distribution for the 
multiplicity of produced particles. However, in order to make a more convincing case for the origin 
of the negative binomial distribution of particle multiplicities, we need a systematic study of bin 
size effects on the extraction of the parameter k from various centralities of heavy ion collisions. 
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